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I.  Introduction 

The  Maximum  Likelihood  (ML)  method  is  widely  regarded  as  the  optimal 
procedure  in  parameter  estimation.  However,  when  applied  to  problems  involving 
composite  signals  and  a  large  number  of  parameters,  the  ML  method  tends  to  be 
computationally  complex  and  time  consuming.  As  an  example,  consider  the 
problem  of  multiple  source  location  estimation  by  signal  observations  using  an 
array  of  spatially  distributed  sensors.  The  received  signals  are  composed  of  the 
contributions  from  the  various  signal  sources,  observed  in  the  presence  of 
additive  noises.  The  ML  estimate  of  the  various  source  location  parameters 
jointly  requires  the  solution  of  a  complicated,  non-linear  optimization  problem 
in  several  unknowns.  Problems  of  that  nature  frequently  arise  in  time  series 
analysis,  array  processing,  radar/sonar,  acoustical  and  geophysical  signal 
processing . 

In  this  report  we  develop  a  computationally  efficient  scheme  for  parameter 
estimation  of  composite  signals,  based  on  the  EM  algorithm.  The  proposed  scheme  is 
optimal  in  the  sense  that  it  converges  iteratively  to  the  exact  ML  estimate  of 
all  the  unknown  parameters  simultaneously. 


II.  Maximum  Likelihood  Estimation  and  the  EM  ATqorithm 

Let  Y  denote  the  data  vector  possessing  the  probability  density 
indexed  by  the  parameter  vector  O  is  a  subset  of  the  Euclidean 

K-space.  Given  an  observed  V,  the  ML  estimate  is  the  value  of  & 

that  maximizes  the  log- likelihood,  that  is 


A/aV 

^^9 


(1) 


Suppose  the  data  vector  Y  can  be  viewed  as  being  incomplete,  and  we  can 
specify  some  "complete"  data  X  related  to  V  by  : 


H  y 


(2) 


where  ^  f .)  is  a  non-invertable  (many-to-one)  transformation.  In  the  multiple 
source  location  problem,  the  'complete'  data  could  be  the  observation  of  the 
various  source  signals  separately,  where  the  observed  (incomplete)  data  is 
the  sum  of  the  signal  contributions  from  the  various  sources. 

The  EM  algorithm  is  directed  at  finding  the  solution  to  (1);  however,  it  does 
so  by  making  an  essential  use  of  the  complete  data  specification.  The  algorithm 
is  basically  an  iterative  method.  It  starts  with  an  initial  guess  f and 
let  be  defined  inductively  by  : 


Max 


(3) 
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where  is  the  probability  density  of  ^  ,  and  f  f»  / 2  r  ^  ' 

denotes  the  conditional  expectation  given  ^  ,  computed  using  the  parameter 
value  ^ The  heuristic  idea  here  is  that  we  would  like  to  choose  0  that 
maximizes  ~  ^  log-likelihood  of  the  complete  data. 

However,  since  Aj available  to  us  (because 
the  complete  data  is  not  available),  we  maximize  instead  its  expectation,  given 
the  observed  data  ^  .  Since  we  have  used  the  current  estimate  rather 

than  the  actual  value  of  0  which  is  unknown,  the  conditional  expectation  is 
not  exact.  Thus,  tne  algorithm  iterates,  using  each  new  parameter  estimate  to 
improve  the  conditional  expectation  on  the  next  iteration  cycle  and  thus  to 
improve  the  next  parameter  estimate. 

The  EM  algorithm  was  first  presented  in  [  11  ,  where  it  is  shown  that 
under  the  usual  regularity  conditions,  the  algorithm  converges  to  the  desired 
result,  that  is  ^  ,  where  each  iteration  increases  the 

likelihood.  The  basic  considerations  leading  to  the  EM  algorithm  (Eq.  (3))  are 
outlined  in  Appendix  A  for  the  convenience  of  the  reader.  The  rate  of 
convergence  of  the  algorithm  is  exponential,  depending  on  the  fraction  of  the 
covariance  of  the  "complete"  data  that  can  be  predicted  using  the  observed 
data.  If  that  fraction  is  small,  the  rate  of  convergence  tends 

to  be  slow,  in  which  case  one  could  use  standard  numerical  methods  to  accelerate 
the  algorithm. 

Ue  note  that  the  EH  algorUhcn  is  not  uhiquely  defined.  The  transformatioh 
H  relating  y  to  y  can  be  any  non-invertable  transformation. 

Obviously,  there  are  many  possible  "complete"  data  specifications  that  will 
generate  the  observed  data.  Thus,  the  EM  algorithm  can  be  implemented  in  many 
possible  ways.  The  final  outcome,  which  is  the  ML  estimate,  is  completely 
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unaffected  by  the  way  in  which  HfO  is  specified  (i.e.,  the  choice  of  "complete" 
data).  However,  the  choice  of  H  {’' )  may  critically  affect  the  complexity  and 
rate  of  convergence  of  the  algorithm,  and  the  unfortunate  choice  of  H  fO  may 
yield  a  completely  useless  algorithm. 

We  shall  proceed  as  follows;  First  we  develop  the  EM  algorithm  for  the 
Linear-Gaussian  case.  This  case  covers  a  wide  range  of  applications.  Then  we 
show  that  for  the  class  of  problems  of  interest  here,  there  is  a  natural  choice  of 
the  "complete"  data,  leading  to  a  surprisingly  simple  algorithm  to  extract  the 
ML  estimates. 
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III.  The  Linear-Gaussian  Case 

Suppose  that  Y'  H  X  ,  where  ^  is  a  matrix  {>n<>i  ),  and  V 

possesses  the  following  multivariate  Gaussian  probability  density: 

•  A 

where  ;^  =  y  if  V  is  real-valued,  A  =  2  if  ^  is  complex  valued, 
and  the  superscript  j¥  denotes  the  conjugate-transpose  operation.  We  shall 
refer  to  this  case  as  the  Linear-Gaussian  case.  Taking  the  logarithm  or  (4) , 
we  obtain: 

■+  ^  (5) 


where  in  C •)  stands  for  the  trace  of  the  bracketed  matrix.  Thus, 

f  f  ''•y  jy  ^  y/-'  y;  ) 

+  ^  j  Wy)/lVf)y'"'-  i  in  f A 'I 't' 


(6) 


where  f ^ j  £  {  if  . 

S  i  nc  e  X  and  Y  are  jointly  Gaussian,  these  conditional  expectations 

are  readily  available  in  the  literature  (e.g.  I  2J  ,  Chap.  5): 

[i- r’^?'">>H7 A 

where  X  is  the  identity  matrix,  and  is  the  “Kalman  Gain" 

defined  by: 

py^) A^?)H''CHA(e'>H*]'^  (9) 

The  EM  algorithm  is  completely  specified  by  Eqs.  (G)  -  (9);  the  algorithm 
iterates  between  calculating  and  ^  5  and  maximizing  the 

expression  in  (6)  with  respect  to  ^  ,  where  each  iteration  increases  the 

likelihood.  We  observe  that  £  { ^0^  ^ y  j 

dependence  on  £  .  Maximizing 

e  (Aj  /x  i 

same  as  maximizing  ^'V  respect  to  £  .  Hence,  the  EM 

algorithm  essentially  requires  the  ML  solution  in  the  X  "lodel  which  might 
be  significantly  simpler  than  the  direct  ML  solution  in  the  ^  model. 
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IV.  Application  to  Signal  Processing 

Let  the  mathematical  model  characterizing  the  observed  signal(s) 
be  given  by 

T.  i  Tf  (10. i) 

where  are  the  vector  signals  received  in  the  absence  of  noise. 

We  stiull  be  concerned  with  composed  of 

K 

('<;«■>  =  5"  (10.2) 

where  1s  the  parameter  vector  associated  with  the  k  signal 

component.  The  specific  problem  we  have  in  mind  is  multiple  source  location 
estimation,  in  which  case  the  Sk  )  are  the  array  signals 

received  from  the  ^  source,  and  are  the  corresponding  source 

location  parameters  (bearing  and  range,  velocity  components,  etc.) 

In  addition,  the  model  given  by  (10)  covers  a  wide  range  of  problems  arising 
in  array  and  signal  processing. 

We  shall  now  consider  the  joint  estimation  of  the  various  9^  for  the 
case  of  deterministic  signals  and  for  the  case  of  stochastic  Gaussian  signals 
separately . 

A.  Deterministic  Signals 

Consider  the  model  of  (10)  under  the  following  assumptions: 

-  The  conditionally  known  to  the  observer  (i.e.  , 

given  9^  ,  we  can  construct  ^ 

h)  ere  vector  zero-mean  white  Gaussian  stochastic  processes 

whose  covariance  matrix  is  E 

where  p  is  a  constant  matrix  and  {EfZ)  is  the  Impulse  function. 
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Under  the  above  assumptions,  the  log-1  llcel  ihood  function  is  given  by 

'  j  j  [  y  ’  1 Hj  if,  )J  ^  fy 

( 

where -c  is  a  normalizing  constant.  The  result  in  (11)  is  a  straightforward 
multi-channel  extension  of  the  known  (deterministic)  signal  In  a  white  Gaussian 
noise  problem  (  [  3) ,  Chap.  4).  In  case  we  are  given  only  a  discrete  set  of 
observations  ^  ^  ^  log-likelihood  is  still  given  by 

(11),  where  the  integral  over  t  is  replaced  by  the  sum  over  the  i;*d  . 


Thus,  the  joint  ML  estimation  of  the  various  obtained  as  the  solution  to: 


In  the  case  of  discrete  observations,  the  corresponding  optimization  problem  is 


I  ;-•/  i-'/  *'■' 

^ 


i/  y 


In  either  case,  the  ML  method  calls  for  the  optimization  of  the  log-likelihood 
function  with  respect  to  all  the  jointly.  Of  course,  brute  force  can 

always  be  used  to  solve  the  problem,  evaluating  the  objective  function  on  a  coarse 
grid  to  roughly  locate  the  global  minimum,  and  then  applying  the  Gauss  method  or 
Newton-Raphson  or  some  other  iterative  gradient-search  algorithm.  However,  when 
applied  to  the  problen  at  hand,  these  methods  tend  to  be  very  complex  and 
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computationdlly  time  consuming. 

Having  the  EM  algorithm  in  mind,  we  would  like  to  simplify  the  optimization 
problem  associated  with  the  direct  ML  approach.  To  apply  the  algorithm  to  the 
problem  at  hand,  the  first  step  is  to  specify  the  "complete*  data.  A  natural 
choice  of  the  "complete"  data  is  obtained  by  decomposing  ^ft)  into 


K 


where 


(14.1) 


(14.2) 

and  the  are  chosen  to  be  mutually  uncorrelated,  zero-mean  Gaussian 

vector  stochastic  processes  satisfying 


Let  the  complete  data  X(t)  be  composed  of  the  various 
from  (14.1),  the  relation  between  the  "complete"  data  and  the  observed 


Then 

( i  nconipl  ete) 


data  is  given  by 

It  )  -  H' 


where 


H-  X] 


J 


(15.2) 


k  terms 
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Sirice  the  "complete"  data  is  Gaussian,  and  the  relation  between  /H)  and 
'i/i)  is  linear,  the  results  developed  for  the  Linear-Gaussian  case  can  be 
applied  here.  The  detailed  derivation  is  given  in  Appendix  B.  The  resulting 
algorithm  is: 

E-step  For  K  ,  compute: 


(16.1) 


M-step  For 


where  the  /5^'d  are  real-valued  positive  scalars  satisfying 

2  •  i  (16.3) 

k  •! 

In  the  discrete  data  case,  the  algorithm  takes  the  form: 

E-step  Compute: 


s'* 


^  /  /  i f  y  "  2  ‘j  ^  ^ 

^•1 


(17.1) 


M-step 


ft 


m:.  I  ft  ■)  -  It  lUj  h 

i‘i 


The  optimization  required  by  (16.2)  is,  in  fact,  the  optimization  problem 
associated  with  the  ML  estimation  of  the  given  separaj^.  observations 

of  the  An  entirely  parallel  statement  can  be  made  with  respect 


to  the  optimization  required  by  (17.2).  Thus,  the  algorithm  decouples  the 
optimization  associated  with  the  direct  ML  approach  into  the  K  separate  ML 
optimizations  as  illustrated  in  Fig.  1 .  Hence  ,  the  complexity  of  the  algorithm  is 
completely  unaffected  by  the  assumed  number  of  signal  components.  As  < 
increases,  we  have  to  increase  the  number  of  ML  processors  In  parallel;  however, 
each  processor  is  maximized  separately.  Since  the  algorithm  is  based  on  the 
EM  method,  it  must  converge  to  the  exact  ML  estimate  of  the  various  9^  'J 
simultaneously,  where  each  iteration  increases  the  likelihood. 

We  note  that  the  ■'d  must  satisfy  the  constraint  (16.3),  but  otherwise 
they  are  arbitrary  free  variables  In  the  algorithm.  The  choice  of  the  does 

not  affect  the  value  of  the  final  estimates;  however,  they  can  be  used  to  control 
the  rate  of  convergence  of  the  algorithm. 

Gaussian  Signals 

Consider  the  model  of  (10)  under  the  following  assumptions: 

-  The  are  mutually  uncorrelated, 

wide  sense  stationary  (WSS),  zero-mean  Gaussian  vector  stochastic  processes 
with  the  spectral  density  matricles  respectively. 

yfH)  is  uncorrelated,  WSS,  zero-mean  and  Gaussian  with  the  spectral 
density  matrix  /V/'tJ)  . 

If  the  observation  interval  T*  -  T;  Is  long  compared  with  the  correlation 
time  (inverse  bandwidth)  of  the  signals  and  the  noises,  the  log-likelihood  is 
given  by  (see  Appendix  C) 

/y  ~  ^  ~ ^  J  ^  ~  ^  j  - ^ ^ ■ 


where 


^  j  C 


(ie.2) 


and 


Tc 


rr 


T. 


k£ 

pC(^j^)=  T  S^r^oj  A/yu))  (18.3) 

/.-/ 

The  summation  in  (18.1)  is  carried  over  all  u)^  in  the  signal  frequency  band. 

In  the  case  of  discrete  observations  ,  the  log-likelihood 

is  still  given  by  (18.1),  where  are  the  discrete 

Fourier  transform  (DFT)  of  .  and  f )  is  given  by 

(18.3),  where  is  the  spectral  matrix  of  the  corresponding  sampled 

signal,  and  /V/’w )  is  the  spectral  matrix  of  the  sampled  noise.  Thus,  in 
either  case,  the  joint  ML  estimation  of  the  various  is  obtained  as  the 

solution  to: 


H 


I  k/i 


I  I  l(o^  de  t  s)  +  )  t )  j 

A  A  y 

Si,---  £, 


A 


(19) 


Thus,  as  in  the  deterministic  signal  case,  the  direct  ML  approach  requires 
the  joint  optimization  of  the  objective  function  with  respect  to  all  the 
Using  the  EM  method,  we  would  like  to  simplify  the  required  optimization. 

Using  the  same  definition  for  the  "complete"  data  as  in  the  deterministic  signal 
case  (Eq.(14))  and  applying  the  results  developed  for  the  Linear-Gaussian  case, 
we  obtain  the  following  algorithm: 
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E-stec 


Compute: 


p'r^,-  i  '■')  j  i/-' )]  + 


#  _-< 


fit (<^, ]  h; ^ » 


M-stei 


it  )r  tn  (B^  g, )  ))] 


(UTl) 


(20.2) 


where  we  define 


(20.3) 


and  the  /3^'d  are  arbitrary  real-valued  scalars  satisfying  the  condition 
(16.3).  The  detailed  considerations  leading  to  Eq.(20)  are  given  in  Appendix  D. 
Perhaps  the  most  striking  feature  of  the  algorithm  is  that  it  decouples  the 
complex  optimization  associated  with  the  direct  ML  approach  into  the  K  separate 
ML  optimizations  as  illustrated  in  Fig.  K  Thus#  as  in  the  deterministic 
signal  case,  we  obtain  a  considerable  simplification  in  estimator  structure  and 
computations. 


Example:  Multiple  Source  Location  Estimation 

The  basic  system  of  Interest  here  consists  of  several  spatially  distributed 
signal  sources  and  an  array  of  several  spatially  distributed  sensors.  The  observed 
signals  can  be  modelled  using  (10),  where  array  signals  received 
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from  the  ^  source,  and  are  the  source  location  parameters  (bearing 

and  range,  velocity  components,  etc.).  Using  the  proposed  scheme  (Fig.  1), 
we  obtain  the  exact  ML  estimates  of  the  various  source  location  parameters, 
while  circumventing  the  complex  multiple  parameter  optimization  associated 
with  the  direct  ML  approach.  For  more  details  we  refer  the  interested  reader 
to  (4  I  and  (  5]  .  To  simplify  the  example,  we  have  considered  the  following 
si tuation: 

-  The  array  consists  of  10  co-linear  and  evenly  spaced  sensors.  The  observed 
data  consists  of  100  time  samples  at  each  sensor  output. 

-  There  are  three  far-field  sources  at  bearings  0.,0.15  and  0.3  radians, 
n\edsured  relative  to  the  Boresight  direction.  Each  source  radiates  a 
triangular  pulse  whose  power  (energy  per  sample)  is  normalized  to  1,  and 
whose  duration  is  20  time  points. 

-  The  additive  noises  are  uncorrelated  from  sensor  to  sensor  with  the  same 
spectral  level  of  0.5  (that  is,  the  post-integration  SNR  per  channel  is 
approximately  16  db). 

We  have  also  assumed  that  LljU'Lt*  ID  ,  where  is  the  sampling  period, 

AL  is  the  spacing  between  adjacent  sensors,  and  U  is  the  velocity  of 

propagation  in  the  medium.  The  problem  is  to  estimate  the  various  source  bearings 
simultaneously. 

In  Fig.  2  we  have  plotted  the  conventional  beam-former  output  as  a  function  of 
bearing.  As  we  can  see,  the  standard  method  cannot  resolve  the  various  sources. 

In  Fig.  3  we  have  plotted  the  outcome  using  our  algorithm.  We  clearly  see  that 
after  about  5  iterations,  the  algorithm  essentially  converges  to  the  true  bearings, 
and  the  three  sources  are  resolved  correctly. 
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V.  Conclusions 

We  have  presented  a  computationally  efficient  scheme  for  parameter  estimation 
of  composite  signals  based  on  the  EM  algorithm.  The  proposed  scheme  is  optimal  in 
the  sense  that  it  converges  iteratively  to  the  exact  Maximum  Likelihood  estimate 
of  all  the  unknown  parameters  simultaneously.  The  method  can  be  applied  to  a 
wide  range  of  problems  arising  in  signal  and  array  processing. 
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Appendix  A  The  EM  Algorithm 

Let  Y  where /Vr*)  is  a  non-invertable  (many-to-one)  transformation. 

Express  densities 

Taking  the  logarithm  on  both  sides  of  (A.l),  we  obtain 

(A.z, 

Denote  y  )and  taking  conditional  expectations  given  Vr  y  for  a 
parameter  value  9*  ,  we  obtain 

fyj  y, .  y  y-- 3 j  t ']  -  urf,  t  v  -  vf!,  t ') 

-  (A. 3) 

Invoking  the  Jensen's  inequality,  we  have  that  Vf  9  ,  Hence 


Eq.(A.4)  forms  the  basis  to  the  EM  algorithm.  The  algorithm  starts  with  an 
initial  guess  ,  and  let  defined  inductively  by 

/W«y  U  (e,  (A.5) 

9 
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Eqs.  (A, 5)  and  (3)  are  identical, 
of  &  tfidt  naxiniizes 
of  the  algorithm  increases  the  value  of 


We  note  that  since  the  value 

then  according  to  (A.4),  each  iteration 
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Appendix  B.  Derivation  of  Eg. (16) 

The  "complete"  data  Xfi)  is  given  by 


y’*) ; 


(B.l) 


where  the  given  by  (14.2),  rewritten  here  for  reference 


(B.2) 


The  (i)  are  chosen  to  be  statistically  Independent,  zero-mean  and 
Gaussian  with  the  covariance  ^ 

where  =  l^l>  ^  ‘  therefore,  mutually  independent  and 

hence,  following  the  usual  considerstlons  ([3],  Chap.  4),  we  obtain 


h  4  -  £  ^03 

^  r  -  f 

hi 

Thus, 


f  f  fy  4v;  ^  )  //-  'jj  e?''^  ’  ]  r 

fc  T5. 

where  contains  all  the  terms  that  are  independent  of  ^  .  The 

are  the  components  of  ^  /tl ,  to  be  computed  from 
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M  XU)  -  ^U)J  f"'  ]  (B.5) 

Using  Eq.(7)  in  (B.4)  and  following  some  straightforward  matrix  manipulations, 
we  obtain  (16.1).  Now,  since  enters  the  right  side  of  (B.4)  only  through  the 

k  term  in  the  sum,  the  joint  maximization  of  (B.3)  with  respect  to  the  various 
decouples  into  the  K  separate  optimizations  as  suggested  by  (16.2). 
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Appendix  C.  Derivation  of  Eg. (18) 

Let  the  data  vector  V  be  denerated  by  Fourier  analyzing  that  is 

f-  )  (c.u 

where  the  are  defined  by  (18.2) .  Since  <dli)  is  Gaussian  and  the 

are  obtained  by  performing  linear  operations  on  'j /i),  then  the 
are  jointly  Gaussian.  Since  both  signals  and  noises  are  assumed  to  be  WSS, 
then  for  observation  interval  T  long  compared  with  the  correlation  time  of 
the  signals  and  the  noises,  the  are  mutually  uncorrelated  (and  thus, 

in  the  Gaussian  case,  statistically  independent)  with  the  covariance  matrix  given  by 

where  is  defined  by  (18,2).  It  follows  that 

Taking  the  logarithm  on  both  sides  of  (C.3)  immediately  yields  (18.1). 


3’*®  Chosen  to  be  statistically  independent,  zero-mean  and 
Gaussian  with  the  spectral  density  matrix  Since  the 

are  assumed  to  be  statistically  independent,  then  the  4  {t  1  are 
statistically  independent,  in  which  case  4nd  statistically 

Independent  whenever  /  yy\  .  Furthermore,  for  long  observation  intervals  and 
WSS  processes,  the  y^/'u^'lat  different  frequencies  are  statistically  uncorrelated 
(and  thus  independent)  with  the  covariance  matrix 


E  f  Xi  (^(  )  (^i  J  5^4  ') 


where  ^ defined  by  (20.3).  It  follows  that 
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{v ''-‘'i?''  -■  ^ IT _ ^  e~  '^1' j  (“-5) 

-  /w  / 

Taking  the  logarithm  on  both  sides  of  (D.5),  we  obtain 


^  #(  -f 

[fo^  Aci  >]  (^.6) 


Taking  the  conditional  expectation  given  Y  s  for  a  parameter  value 


we  obtain 


l^V'i 


where  is  the  block  matrix  of  computed  from 

Y )  --  5 / //'«, )  yL^ ) /HYn ' --  y<^t  );i''”} 


(D.8) 


Using  Eq.(8)  in  (0.8)  and  following  some  straightforward  matrix 
manimpulations ,  we  obtain  (20.1).  Since  the  a  C  )  term  of  the  double  sum  on 
the  right  side  of  (0.7)  depends  only  on  ,  then  the  joint  maximization  of 

(0.7)  with  respect  to  the  various  decouples  into  the  1^  separate 
optimizations  as  suggested  by  (20.2). 
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Flgure  Captions 

Fig.  1:  Signal  Processing  Via  the  EM  Algorithm. 

Fig.  2:  Conventional  Beamforming. 

Fig.  3;  Multiple  Source  Location  Estimation  Via  the  Proposed  Algorithm. 
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Fig.  3:  Multiple  Source  Location  Estimation  Via  the  Proposed  Algorithm. 
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